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Abstract 

^ _ Wave propagation in a stratified fluid / porous medium is studied here using 

analytical and numerical methods. The semi-analytical method is based on 
■ an exact stiffness matrix method coupled with a matrix conditioning pro- 

cedure, preventing the occurrence of poorly conditioned numerical systems. 
Special attention is paid to calculating the Fourier integrals. The numeri- 
cal method is based on a high order flnite-difference time-domain scheme. 
Mesh reflnement is applied near the interfaces to discretize the slow com- 
pressional diffusive wave predicted by Blot's theory. Lastly, an immersed 
interface method is used to discretize the boundary conditions. The numeri- 
cal benchmarks are based on realistic soil parameters and on various degrees 
of hydraulic contact at the fluid / porous boundary. The time evolution of 
the acoustic pressure and the porous velocity is plotted in the case of one and 
four interfaces. The excellent level of agreement found to exist between the 
two approaches conflrms the validity of both methods, which cross-checks 
them and provides useful tools for future researches. 
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1. Introduction 

This study focuses on the propagation of transient 2D mechanical waves 
emitted by a source point in a fluid half-space over a stratified poroelastic 
medium. This configuration occurs in situations such as those encountered 
in the field of underwater acoustics and civil engineering. Natural or arti- 
ficial media, presenting unidirectional varying properties, can be modelled 
as multilayered structures. The wave propagation is described by means 
of acoustic equations in the fluid and low-frequency Biot's equations in the 
poroelastic media Hydraulic exchanges between the fluid and the first 
porous layer are modeled using various boundary conditions [2I, [sf. Classical 
textbooks such as 0, Isf can be consulted for a detailed analysis of Biot's 
equations (which involve a fast compressional wave and a shear wave, as in 
elastic media, and a slow compressional wave). In the case of a single fluid 
/ porous interface, many theoretical studies have dealt with the reflection 
/ transmission coefficients [sj and with the properties of the surface waves 
0, S M, EdI 11 1 . The aim of the present study is to solve this problem in 



the case of an arbitrary number of layers, using two radically different ap- 
proaches: a semi-analytical approach and a purely numerical one. 

On the one hand, various analytical methods have been developed in the 



case of a single fluid / porous interface: Cagniard-de Hoop's method [l2|, [13 



for an inviscid medium, and [IJ, |15| for a viscous saturating fluid, to cite but 
a few. In situations involving a larger number of interfaces, the strategies 
usually adopted are based on transfer matrix, stiffness matrix, or transmis- 
sion and reflection matrix methods. These approaches were first developed 



for use with viscoelastic media [16|, ll7|, ll8| , and only a few studies using mul 



tilayer approaches to poroelastic models have been published so far because 
of the complexity of poroelastic models and the poor conditioning of the re- 
maining matrices. The exact stiffness matrix method has been applied in the 



context of poroelasticity to harmonic 2D problems [19|], transient 2D prob 



lems (20j, and problems involving axisymmetric quasi-static configurations 



21 . All the latter authors stressed the numerical difficulties encountered 



due to poorly conditioned matrices. Methods have been proposed in [2^ for 
eliminating the mismatched terms using the transmission and reflexion ma- 
trix approach in the case of axisymmetric geometries, and recently, using the 



2 



exact stiffness matrix method in 23|. This method is extended here to the 
couphng between ffuid and porous media by simulating the various interface 
conditions pertaining at the boundary between the fluid and the stratified 
porous media. 

On the other hand, only a few numerical methods have addressed realistic 
acoustic / poroelastic configurations: for this purpose, a spectral-element 
method 2J], a pseudospectral method 25| and a discontinuous Galerkin 
method 26f have been developed. Three difficulties have to be overcome 
to be able to perform reliable and efficient numerical simulations. First, 
the slow compressional Biot wave is highly attenuated, which drastically 
affects the stability of any explicit numerical scheme 27|. Secondly, the slow 



compressional wave is always located near the interfaces, and thus plays a 
key role in the balance equations; discretizing this small-scale wave therefore 
requires a particularly fine spatial mesh. Thirdly, large numerical errors are 
usually introduced by the boundaries due to the poor discretization of their 
geometrical and physical properties, and also due to issues involved in the 
numerical analysis. In the present study, it is proposed to adopt a strategy we 



have described in previous articles (28|, |29| . A fourth-order finite-difference 



scheme with splitting is used to integrate the evolution equations, so as to 
optimize the integration time step. Specific solvers are used in the case of the 
fluid and the porous media. Their coupling is obtained using an immersed 
interface method to discretize the boundary conditions, which gives a subcell 
resolution of the geometries. Lastly, a space-time mesh refinement procedure 
is applied around the interfaces to account for the small scale of the slow 
compressional wave. 

The paper is organized as follows. In section [2l notations, governing equa- 
tions and boundary conditions are stated. In section [31 the semi-analytical 
approach is described, focusing on the coupling between fluid and porous 
media and on the integration of oscillating functions. In section HJ the key 
components of the numerical approach are recalled. In section |5l numerical 
experiments are performed on situations involving one and four interfaces, 
with realistic soil parameters and under the low-frequency Biot theory. Com- 
parisons between the analytical and numerical solutions obtained confirm the 
validity of both approaches. From the numerical point of view, this study 
provides a useful benchmark for future studies on viscous saturating fluids 
that have not been available so far 2J, |25|, |26|. From the semi-analytical 



point of view, these comparisons also show the relevance of using the chosen 
integration parameters, for which no theoretical criterion exists. Concluding 
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remarks and some future lines of research are presented in the last section Ei 

2. Statement of the problem 

2.1. Governing equations 
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Figure 1: Configuration under study. Fluid medium over a multilayered porous medium. 
Source point ys in the fluid. 

The 2D configuration under investigation is a fluid half-space VLq over a 
stack of homogeneous poroelastic layers VLj (j = 1, ■ ■ ■ , A^), as shown in Fig. 
[U The X and y geometrical axes point rightward and upward, respectively. 
The plane and parallel interfaces are located at ?/ = ctj < 0, where = 0. 
A source point at {xg = 0, ys > 0) in the fluid emits cylindrical waves. 

In the fluid domain Qq, the physical parameters are the density pj and 
the celerity of acoustic waves c. The acoustic equations are written as follows 

p = -pc^ V.U, 

1 (1) 

Ap- -^p= -S{t)6{x) 5{y - ys), 

where U = {Ux, UyY is the fluid displacement, p is the acoustic pressure, 
and S{t) is a causal source term. The overlying dot denotes the derivative in 
terms of time t. 
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The poroelastic media flj (j = 1,---,A^) are modeled using the low- 
frequency Biot theory [H, 0]- The physical parameters are 

• the dynamic viscosity r] and the density p/ of the saturating fluid. 
The latter is assumed to be the same as in Qq, and the notation p/ is 
therefore used in both cases; 

• the density ps and the shear modulus fi of the elastic skeleton; 

• the connected porosity (p, the tortuosity Coo, the absolute permeability 
K, the Lame coefficient of the dry matrix Aq, and the two Biot coeffi- 
cients P and m of the isotropic matrix. 

We thus obtain the total density p = 0p/ + (1 — </>) p^. The low-frequency 
Biot model is valid at frequencies below the critical value defined as 

27raoo i^Pf 

Based on the constitutive equations and the conservation of momentum in 
porous media, one obtains |l| 

a = (AoV.u-/3p) 1 + 2 fie, 
p = —m (/3 V.u + V.w) , 

(3) 

V a = pu + pj w, 

— Vp = pf u H ■ — w H — w, 

where u = (u^., UyY is the solid displacement, U = (^4, UyY is the fluid 
displacement, w = (U — u) = {w^, WyY is the relative displacement, I is 
the identity tensor, a is the stress tensor, e = | (V u + V* u) is the strain 
tensor, and p is the pore pressure. 

2.2. Boundary conditions 

The governing equations ([T]) and have to be completed by a set of 
boundary conditions along the N plane interfaces y = aj (j = 0, ■ ■ ■ , — 1). 
For this purpose, we take [g]j to denote the jump in a function g from Qj to 
flj+i across aj 



The fluid / porous interface is modeled with the following boundary con- 
ditions 0, B H 0, H 



yJo 



[w. 



yyjo 



[P]0 = ^ {Wy)o 



(5a) 
(5b) 
(5c) 

(5d) 



Eq. fl5dp involves the hydraulic permeability of the interface /C, resulting in 
the following limit cases: 

• if /C — 7- +00, then Eq. (l5d|) becomes [p]o = 0, describing open pores; 



if /C — 0, then no fluid exchange occurs across ao, and Eq. fl5dp is 
replaced by {'Wy)Q = 0, describing sealed pores; 

if < /C < +00, then an intermediate state between open and sealed 
pores is reached, describing imperfect pores. 



The porous / porous interfaces aj {j 
perfect bonded contact 0,01: 



1, ■ ■ ■ — 1) are assumed to be in 
njj = 0, [uy]j = 0, [Wy]j = 0, [a^^y]j = 0, [ayy]j = 0, [p]^ = 0. (6) 



3. Semi-analytical solution 

3.1. Helmholtz decompositions 

Pressure and stress components are eliminated from Eqs. ([3]), giving a 



[u, w) second-order wave formulation 

V7^V7 „\ , ,. V72 



(Ao + ;U-Fm/3^) V(V.u) + /i u + m /3 V(V.w) = pii + p/w, 
m/3 V(V.u) +mV(V.w) =p^u + + - w. 

K 



(7) 



The solid and relative displacements in Eqs. can be expressed by 

u = V<^ + Vx*, w = Vv9'' + Vx*^ (8) 



6 



where if and (f"^ are scalar potentials, and ^ and are vector potentials. 
One introduces mass, stiffness and damping matrices 



(9) 



\q + 2n + m (3'^ m (3 \ //uO 

Kp = I ^ , Ks = 

m p ml \ 

P Pf \ / 

— / V 

The Fourier transforms in time t and space x are defined by 

/H-cxD -1 r-\-oo 

rico) exp{-tcot)dco, riu) = — / fit) exp{tcot)dt, 

/-\-00 2 /' + 00 

g{k^) exp{ik^x)dk^, g(k:,) = — g{x) exp{-ik^x)d. 
-co J— oo 

(10) 

Applying these x and t Fourier transforms to relations ([7j) and ([8]) yields de- 
coupled ordinary differential systems in the frequency-wavenumber domain, 
associated with fast and slow compressional waves Pf and Pg, and with shear 
wave S, respectively 



X. 



(11) 



Taking the Fourier transforms of Eqs. ( uTj) over y gives the dispersion rela- 
tions 

det (4^- Kp - w^M - z w C) = 0, kl = —{p + pfG), (12) 

A'" 

where kpj [j = /, s) are the wavenumbers of the fast (/) and slow (s) 
compressional waves, and ks is the wavenumber of the shear wave, see 
Appendix A, Eq. (lA.ip . G is defined in Appendix A Eq. flA.2|) . This 



gives the phase velocities cpj = u / ^{kpj) and cs = uj/^{ks), where 
cpf > max(cp^, Cs). 
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3.2. Exact stiffness matrix approach 

The exact stiffness matrix approach is based on vectors of transformed 
displacement and stress components 



20 , defined as 



^Uy, IWy) , 



yy 



Analytical expression for the transformed displacement is then solution of 
the matrix system 



,Q^Z(/z,) 





(6x6) 



where 



a 



■ a. 



is the thickness of layer 17^ (j = 1, ■ • • , ■ 



(13) 

1). Further 



information about the matrices 'Zi{hj), S^'^, Q"^'^ is given in Appendix A 



Eqs. flA.Sp . flA.4p and f lA.Sp . respectively. The superscripts T and R stand 
for transmitted and reflected waves respectively. Special attention is paid to 
the conditioning of the matrices. Increasing exponential terms corresponding 
to reflected waves are excluded from the formulation. These are included in 
the unknown wave potential vectors, which do not need to be calculated, 



see [23|. Only decreasing exponential terms are therefore explicitly included 
in the expressions of Z{hj), Appendix A Eq. (lA.Sp . Based on Eqs. ([6]), 
a classical assembling procedure between the porous layers is used, which 
involves the continuity of the stresses and displacements at each interface. 
In the case of the half-space [hs) Qn under the layers, only the transmitted 
terms are kept and T'^^ = S^(Q"^)~^. The whole matrix system involving all 
the poroelastic subdomains has dimension 3 A^ x 3 A^. 

3.3. Assembling of the fluid / porous interface 

In the fluid domain Qq, the acoustic pressure p* is the sum of an incident 
cylindrical wave emitted by the source point and a reflected wave 



— * 

p 



p + p 



iS{u) exp{ikyf \ y - ys \) 



4n 



k 



yf 



+ A exp(i kyfy), 



(14) 



where A is determined by the boundary conditions 1^ along oq, and kj = 
uj'^/cl = kl + k'^r. Taking kyj with '^{kyf) ^ 0, the pressure vanishes at 



y — 7- +00 (the radiation condition) 



yf) _ _ 

Expressions for and Uy are easily 
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deduced from Eqs. 0141) and ([T]). In the matrix condensed form, the part to 
be included into the overall system is written as follows 



/ ipfUJ 



I PfLO 



V 



i pf uj^ 



yf 



I PfUJ 



yf 



[ \ 

Ic ) 



I [u 



I [w 



( . ^_,^+ , S{ul expjt kyf Vs) \ 



Coupling system (fT3|) with (fT5|) yields the overall system 



•A 



-^22 



-^32 



I PfUJ 

i PfUj"^ 



-^23 



I PfUJ 



''yf 



-^33 



Pf^ 

Kf 



Ic 



\... 



/(K)o\ 



I [W 



2tx kyf 
_ Sju) exp{ikyf Us) 

2TX kyf 

(15) 

/ \ 

S{uj) expi^ikyfi/s) 



.../ 



V ••• / 



27r 



k. 



yf 



S{uj) exp{ikyfys 



271 



^yf 



\ 



J 



(16) 

The solution to system (fT6|) gives the horizontal, vertical solid and relative 
displacements at each interface. The displacements, stresses, velocities and 
acoustic pressure in each domain VLj can then be obtained. 

3.4- Processing the oscillating integrals 

Carrying out inverse Fourier transforms over the horizontal wavenumber 
kx and the angular frequency uj requires the use of numerical integration 
procedures. 

Two problems have to be overcome when calculating the integral over the 
horizontal wavenumber. First, the integrand shows fast oscillatory behavior 
due to the factor exjg{ikxx). Secondly, the envelope of the maximum am- 
plitudes shows sharp peaks, as shown in Fig. [2l Alternative techniques to 
the usual Fast Fourier Transform for evaluating the integral properly can be 



found in [3l|, |32|, |33| in the context of elastodynamics, and |19l . |20| in that of 



poroelastodynamics, for example. Here we use the adaptive Filon quadrature, 
which is particularly accurate and suitable for dealing with these integrals. 
The transformed fluid pressure, pore pressure, stress tensor normal compo- 
nent, vertical displacement and velocity are symmetric terms with respect to 
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120 



Figure 2: Processing the integrals. Real part (left) and imaginary part (right) of 
iiy (kx, t^o) (blue) and Uy (k^, wo) cos(A:a; a;) (red). This signal is obtained at recorder 
i?4 in the case of = 4 interfaces, with sealed pore conditions; see section [5] for details. 



kx- These functions, which are denoted by h , can be written as follows 
h{x,y,t)= / h {kx, y, oj) cos{kxx) dkx? exp{—iut) du. (17) 



In line with the Filon quadrature theory, the integrand is separated into 
the product of a slowly-varying component and a fast-oscillating one. The 
oscillatory kernel is rigorously accounted for, since only the complex function 
h {kx, y, uj) is replaced by a (second degree here) polynomial approximation. 



The resulting integral is then computed analytically [34 . 

One possible adaptive procedure consists in dividing the entire interval 
into several parts based on what is known about h {kx, y, ou). Because of 
the sharp changes in the integrand occuring around the wavenumbers of the 
propagating waves, the wavenumbers are calculated and sorted out. The 
quadrature is performed by discretizing finely in the neighborhood of these 
wavenumbers and more coarsely farther away. The integral is truncated 
depending on the highest wavenumber and adapted to each frequency. 

Concerning the numerical processing over u, the physical space-time do- 
main h{x, y, t) term is real. We therefore obtain 

h{x, y, t) = I <^2 I h {kx, y, u) cos{kxx)dkx > exp{—iut) j du. 

(18) 



\ \^ JO 
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Numerical integration is done using a Simpson quadrature, and the Nyquist- 
Shannon sampling theorem is checked. 



4. Numerical solution 



4.I. Numerical scheme 

Velocity-stress formulation for the governing equations can be obtained 
from Eqs. ([T]) and Setting 



A 



Ux, Uy, P 



(Uxi My, Wxi '^yi ^xxi ^xy) '^yyi 



in Qq 
in Q 



'J' 



J 



(19) 

gives the first-order linear system of partial differential equations with source 
term 

A + A — A + B — A = -CA + F. (20) 

ox ay 

In Eq. (120|) . A, B and C are 3x3 matrices in Qq, and 8x8 matrices in Qj 
{j = 1, - ■ ■ ,N); the vector F accounts for the acoustic source in Eqs. ([1]). In 
Qq, C = 0. Due to the non-zero diffusive term C present in the poroelastic 
media flj, an unsplit discretization scheme for Eq. (I2U|) would not be suitable, 
leading to a penalizing time step. A much more efficient method, which was 
used in 28|, consists in using a Strang splitting procedure. The following 
partial differential equations are solved successively 



d d 
A + A — A + B — A 

ox ay 



0, 



A 



A 



CA, 



(21a) 
(21b) 
(21c) 
on a 

uniform Cartesian grid, with spatial mesh sizes Ax, Ay and a time step At. 
This explicit two-step finite-difference scheme is accurate to the fourth-order 
in both time and space, and satisfies the optimum CFL stability condition. 
The diffusive part (I21bp is solved exactly in the poroelastic layers fij, see Eq. 



The propagative part (12 lap is solved using an ADER 4 scheme 35 



m 



28 



This step induces no stability restrictions. The fluid VLq has no 
diffusive part, because there is no attenuation. Lastly, numerical integration 
is performed on VLq to account for the source term (j21cp . 
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4-2. Mesh refinement 

In the low-frequency range, the slow compressional wave is a diffusive 
non-propagating solution with very small wavelength Xp^. A very fine spatial 
mesh is therefore required. Since this wave does not propagate and is always 
located near the interfaces, space-time mesh refinement is a good strategy. 



For this purpose, a steady-state version of the algorithm presented in [36 
is developed. Each interface is inserted into a refined grid, where both the 
spatial meshes and the time step are divided by the refinement factor q. 
This procedure ensures that the same CFL number will be obtained in each 
grid. The coupling between coarse and fine meshes is achieved by performing 
spatial and temporal interpolations. The factor q is used to obtain the same 
discretization of the slow wave on the refined zone, as with the fast wave 
on the coarse grid. Taking /o to denote the central frequency of the signal, 
direct calculations give g(/o) = cp/(/o)/cp^(/o). 

4-3. Immersed interface method 

The discretization of the boundary conditions requires special care, for 
three reasons. First, if the interfaces do not coincide with the uniform mesh- 
ing, then geometrical errors will occur. Secondly, the conditions (E]) and 
(E]) will not be enforced numerically by the finite-difference scheme, and the 
numerical solution will therefore not converge towards the exact solution. 
Lastly, the smoothness of the solution required to solve Eq. fl20|) will not be 
satisfied across the interface, which will decrease the convergence rate of the 
ADER scheme. 

To overcome these drawbacks without detracting the efficiency of Carte- 



sian grid methods, an immersed interface method [28|, |29| is used. The latter 
studies can be consulted for a detailed description of this method. The basic 
principle is as follows: at the irregular nodes where the ADER scheme crosses 
an interface, modified values of the solution are used on the other side of the 
interface instead of the usual numerical values. 

Calculating these modified values is a complex task involving hi gh-o rder 



derivation of boundary conditions (l5])-(|6]), Beltrami- Michell relations |37[ and 
singular value decompositions. Fortunately, all these time consuming proce- 
dures can be carried out during a preprocessing stage and only small matrix- 
vector multiplications need to be performed during the simulation. After 
optimizing the code, the extra CPU cost can be practically negligible, i.e. 
lower than 1% of that required by the time-marching procedure. In addi- 
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tion, by choosing the order of the interface method carefully, it is possible to 
achieve the overall fourth-order accuracy of the ADER 4 scheme. 



5. Numerical results and discussion 

5.1. Configuration and physical parameters 



Saturating fluid p/ (kg/m^) 
c (m/s) 
Tj (Pa.s) 

Grain ps (kg/m^) 

/i (Pa) 

Matrix (p 

K (m^) 
Ao (Pa) 
m (Pa) 



Phase velocities 



Cp/(/o) (m/s) 
cpsifa) (m/s) 
csUo) (m/s) 
est (m/s) 
fc (kHz) 



1000 

1500 

10-3 

2760 

3.40 10^ 

0.24 

2.3 

3.910^^3 

3.5 10^ 
8.1 10^ 
0.855 



2636.8 

571.5 

1210.8 

1007.1 

42.5 



2i 



1000 
1500 
10-3 
2644 

7.04 10^ 
0.20 
2.4 

3.610-^3 

5.5 10^ 
9.7 10^ 
0.720 



3263.7 

649.5 

1751.0 

36.8 



Table 1: Poroelastic media 0,2 j-i and 1^2 j {j 
acoustic properties at /o = 20 kHz. 



1, • • • , N): physical parameters and 



Two configurations are proposed for testing the above semi-analytical and 
numerical methods. The influence of the various fluid / porous boundary 
conditions on the computed fields is also investigated. The first test involves 
a single interface (A^ = 1) between a fluid and a porous medium. The second 
test is performed on a layered medium involving 4 interfaces (A^ = 4). The 
accuracy and robustness of the semi-analytical and numerical approaches are 
proven in the case of realistic soil parameters. 

The acoustic medium VLq is water (p/ = 1000 kg/m^, c = 1500 m/s). The 
poroelastic layers consist alternately of water-saturated sand in VL2j-i jsj, 
and Berea sandstone in VL2j 23|. The poroelastic properties of these media 
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Figure 3: Time evolution (left) and spectrum (right) of the source S defined by Eqs. p2|l 
and (|23p . The vertical dashed line gives the central angular frequency loq = 2t: /q. 

are summarized in Table [TJ In the case where an imperfect hydrauhc contact 
occurs along the interface ao, the hydraulic permeability is /C = 5. 10~® 
m/s/Pa. 

The time evolution of the source S appearing in Eq. ([T]) is a combination 
of truncated sinusoids 



where Pm = and = 27r /o; the coefficients a™ ensuring smoothness 
of the solution are: a\ = 1, a2 = —21/32, 03 = 63/768, 04 = —1/512. The 
Fourier transform of Eq. f l^2]) is 



The central frequency in Eqs. f l22|) and fl23|) is /o = 20 kHz, which is much 
smaller than the critical Biot frequency of poroelastic layers, see Eq. ([2]) and 
Table [TJ The time evolution and Fourier transform of S are presented in Fig. 
El In the following test cases, the source is always applied in the fluid domain 
at point [xs = 0, i/s = 4.10"^) m. 

5.2. Test 1: N = 1 interface 

The interface between water and saturated sand is set at ao = m. With 
the present numerical method, the whole computational domain [—2,2]^ m^ 




/o 



1 



(22) 




(23) 
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x(m) 

Figure 4: = 1 interface with open pores. Snapshot of p (in the fluid) and a^x (in the 
poroelastic medium) at t = 0.72 ms. Green-red: compressional waves; yellow-magenta: 
shear wave. 
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is discretized using a regular Cartesian grid with Ax = Ay = 2.10"^ m. This 
grid is refined locally by a factor g = 5 in the neighborhood of the interface in 
order to satisfy the criterion defined in section 14.21 With these parameters, 
the fast and slow compressional waves are discretized by at least 75 grid 
points per wavelength, and acoustic and shear waves by about 37 points per 
wavelength. 

The numerical solution obtained {p in the fiuid, a^x in the poroelastic 
medium) is shown in Fig. HI at time t = 0.72 ms under open pore conditions. 
Since the source is located very near the interface, the refiected acoustic wave 
in Qq cannot be distinguished from the circular incident wave; the headwaves 
are also clearly visible. The transmitted fast compressional wave (circular 
wavefront) and shear wave (more complex structure) can be seen in Qi. 
Pseudo-Stoneley waves can also be observed along the interface, just behind 
the shear wave. The transmitted slow compressional wave cannot be seen 
here; it remains located along the interface on very small spatial scales. 

In this configuration, the semi-analytical results were obtained with the 
following set of numerical parameters. As previously explained in section 
13. 4^ wavenumbers corresponding to the propagating waves are calculated and 
sorted out. Let kx,min and fcx.max denote the lowest and highest wavenumbers 
of propagating waves. The integral over kx is divided into three parts. The 
second part includes the sharp evolution of the integrand, and integration is 
performed in this part with a finer grid. We perform 

1. integration over G [0, 0.75 A;a;^min] with 500 subintervals; 

2. integration over kx G [0.75 fca^^min, 1-25 /c^^^max] with 5000 subintervals; 

3. integration over kx G [1.25 fc^^^max; 20 /c^^^max] with 200 subintervals. 

The numerical integration over u is led with 2000 regular subintervals, up 
to the maximum angular frequency Wmax = 710^ rad.s"^. Sampling with 
respect to time was performed with 4000 points, where t G [0, 0.70] ms in 
Fig. El 

The time history of the pressure p and of the vertical velocity Uy obtained 
with these two methods is presented in Fig. El The receivers were placed 
near the interface in order to capture the surface waves. The pressure was 
recorded in the fiuid domain at point Rl {x = 0.5, y = 0.02) m while the ver- 
tical velocity was recorded in the porous medium at R2 [x = 0.5, y = —0.012) 
m. Results were normalized with respect to the maximum value of the pres- 
sure (Rl) or velocity (R2) respectively. Excellent agreement was obtained 
between the semi-analytical and numerical values, under all the interface 
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Figure 5: = 1 interface, with open pores (top), imperfect pores (middle) and scaled pores 
(bottom). Time evolution of the acoustic pressure p at recorder Rl {x — 0.5, y = 0.02) m 
(left) and that of the vertical velocity iiy at recorder R2 {x = 0.5, y = —0.012) m (right). 
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pore conditions tested. The validity of the numerical strategy used to com- 
pute the integrals in the semi-analytical method was therefore confirmed. 
From the physical point of view, Figure |5] presents three peaks due to the 
following waves: Pf for the fast compressional wave {tpf ^ 0.21 ms), F for 
the direct fluid wave (tp = 0.33 ms), and St for the pseudo-Stoneley wave 
(tgj ^ 0.51 ms). The influence of the pore conditions adopted at the fluid 
/ porous interface can be clearly seen in Fig. [61 Note that the direct fluid 
wave occurring in the acoustic pressure component (Rl) was not affected by 
the pore conditions, contrary to what occurred with all the other peaks. 

5.3. Test 2: N = 4 interfaces 

The computational domain was the same here as in test 1 and the co- 
ordinates of the 4 interfaces were set as follows: = m, ai = —0.1 m, 
a2 = —0.15 m, and = —0.6 m, with layer thicknesses ranging from 0.5 to 4 
wavelengths of the fast compressional wave. Once again. Ax = Ay = 2.10"'^ 
m were used for the coarse grid, and Ax/5, Ay/ 5 in the more highly refined 
grids surrounding each interface. 

The numerical solution is presented in Fig. [TJ where complex interactions 
and wave conversions can be observed. The time evolution of the acoustic 
pressure at point R3 (x = 1, ?/ = 0.02) m and that of the vertical velocity Uy 
at point R4 [x = 1, y = —0.112) m are shown in Fig. [HI Once again, excellent 
agreement was observed between the semi-analytical and numerical values 
obtained. With the conditioning matrix technique, no oscillation occurs in 
the semi-analytical values, and the layers therefore do not need to be divided 
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Figure 7: = 4 interfaces involving open pores at the fluid/porous boundary. Snapshot 
of p (in the fluid) and a^x (in the poroelastic layers) at t ~ 0.69 ms. 
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into sublayers, which makes the calculations more efficient 23 



Differences between the pressures computed under open, imperfect and 
sealed conditions can be clearly seen in Fig. |9l especially in the successive 
diffiacted waves. This shows the need for an accurate means of treating the 
pore conditions whatever the numerical method used, even in the case of a 
single fluid / porous interface. 



6. Conclusion 

The propagation of transient waves in a stratified 2D fluid / poroelastic 
medium was studied here using two different approaches: a semi-analytical 
and a numerical approach. On the one hand, the analytical solution is fast. 



and can therefore be efficiently plugged into imaging codes [38|, |39[. This 
approach can only be used with academic geometries, however. On the other 
hand, the numerical solution can be used to handle much more complex 
geometries and media, but the full space-time grid is not suitable when the 
solution is sought only at some receivers. 

Both methods are also able to deal with various boundary conditions, 
which can have non-negligible effects on the various fields, as we have seen 
here. As far as we know, the case of imperfect pores has never been treated 
previously using a numerical approach. Since the assumptions underlying 
the present semi- analytical and numerical approaches are radically different, 
the comparisons between the results obtained constitute an authentic cross- 
check, and the results obtained therefore provide useful benchmarks for future 
researches. 

We are currently extending this study to include frequencies greater than 



Eq. ([2]), using the JKD model |40|. The frequency correction can be straight- 
forwardly incorporated into the analytical method, but the time-domain 
numerical modeling needs to be considerably adapted because of fractional 



derivatives involved |4l| . Work is currently under way on these lines [42 
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Figure 8: = 4 interfaces, with open pores (top), imperfect pores (middle) and sealed 
pores (bottom). Time evolution of the acoustic pressure p at recorder R3 {x = I, y = 0.02) 
m (left) and that of the vertical velocity iiy at recorder R4 (x = 1, y = —0.112) m (right). 
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Figure 9: = 4 interfaces. Comparison between the different pore conditions: pressure 
(left) and vertical velocity (right). 



Appendix A. Notations used along section [3] 

The vertical components of the wavenumbers (|T2l) . which are denoted by 
an y index, satisfy 

^y-Pj ~ ^"pj ~ ^yS — ~ ^x- (A-l) 

To ensure the radiation condition, we choose Q{kyPj, kys) > 0. Then, the 
coefficients associated with compressional waves and the shear wave are 



{kyp,: + k^)/3m- uj^ Pf 

Fpj = —7j2 -TT\ \ — 2 TT^- J = /,s (fast, slow), 

-{kypj + k;)m + u^aoo Pf/(l) + iUJr]/f^ 



ttoo Pf KUJ + if] (j) 

(A.2) 

The matrices Z, S'^'"'^ and Q'^'-'^ in system ( !T3|) are given by 

Z{hk) = diag < exp(? kypf hk), exp{i kyPs hk), exp{i kys hk) >, (A.3) 
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where diag stands for the diagonal matrix, 
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Sji — 2 fi kypf kx, 
= 2 /i kyPs kx, 
^13 = {kyg — kl), 
Sli = % {-{kip J + (Ao + m /32 + m 

q2 



0I2 = ^ {-{k^Ps + kl) (Ao + + m 

*^23 = 2 2 kyS kx , 

Sl = -zm{klp^ + kl) {Fpj + P) 
Sg = -tm{k'+kl) {Fps + P), 
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i k. 



ys 



[ Qli 



\ kyPfFpf kyPsFps —kxG j 



21 



\ -Q 



31 



Q12 
~Q^2 
~Q32 
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